% save('data/experiment-60min.mat', 'inst_count_vec', 'rt_vec');
% load('data/experiment-60min.mat');

%  clear all
%  clear
%  load('data/experiment-60min.mat');
%  open policy/test_fake_lqr.m
%  test_fake_lqr().build_model( inst_count_vec, rt_vec)

            %----------
            n_st=1;% 1 state variable
            n_in= size(inst_count_vec,2);  %11 inputs
            n_out= 1; %size(rt_vec,2); % 1 output
            n_env=1; % environmental disturbances
            A=ones(n_st,n_st);
            B=ones(n_st,n_in);
            K=ones(n_st,n_env);
            C=eye(n_out,n_st);
            D=zeros(n_out,n_in);
            
            m_rt_vec=0; %2
            [  d_rt_vec]= obj.linearize_around(rt_vec , m_rt_vec);  %  2  seconds desired
            m_inst_count_vec=0; %2
            [ d_inst_count_vec] = obj.linearize_around(inst_count_vec,...
               ones(1, size(inst_count_vec,2))* m_inst_count_vec);
            
            u =  d_inst_count_vec;
            y = d_rt_vec;
            
            data=iddata(y, u , obj.samplingTime);
            % m=idss(A,B,C,D,K,'Ts', obj.samplingTime,'DisturbanceModel','Estimate','InitialState','Estimate');
            m=idss(A,B,C,D,K,'Ts', obj.samplingTime,'DisturbanceModel','None','InitialState','Estimate');
%             m0.As = [NaN,0;0,NaN];
%             m0.Bs = [NaN;NaN];

            m.Cs = C;
            m1=init(m);
            mm=pem(data,m1);
            
            % Design LQ-optimal gain K            
%         K = lqry(sys,10,1)  % u = -Kx minimizes J(u)
            % Add integrator state dz/dt = -phi
            A_aug = [0 -1;
                            0 sys.a];
            B_aug = [zeros(1,11) ; sys.b];

            % LQR gain synthesis
            Q = blkdiag(1, 0.1); % for integrated difference and for crappy state (here response time) 
            R = diag(prices_amz().get_hourly_price);
            lqr_gain = lqr(A_aug,B_aug,Q,R);
            
            